Phases of a rotating Bose-Einstein condensate with anharmonic confinement 
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We examine an effectively repulsive Bose-Einstein condensate of atoms that rotates in a quadratic- 
plus-quartic potential. With use of a variational method we identify the three possible phases of 
the system (multiple quantization, single quantization, and a mixed phase) as a function of the 
rotational frequency of the gas and of the coupling constant. The derived phase diagram is shown 
to be universal and the continuous transitions to be exact in the limit of weak coupling and small 
anharmonicity. The variational results are found to be consistent with numerical solutions of the 
Gross-Pitaevskii equation. 



PACS numbers: 03.75.Hh, 03.75.Kk, 67.40.Vs 
I. INTRODUCTION 

When rotated, superfluids are expected to show a 
structure of quantized vortex states. This effect has 
been confirmed experimentally in Bose-Einstein conden- 
sates of atoms which typically are confined in harmonic 
potentials [f-5]. In contrast to traditional experiments 
on superfluid liquid helium contained in a "bucket" , the 
presence of a harmonic potential introduces the trap fre- 
quency, CO, which sets the scale for the critical frequency 
of rotation of the gas fie for vortex formation. The same 
frequency w also sets an upper bound for the rotational 
frequency of the atoms. When Q = lu, the centrifugal 
force exactly cancels the confining force, and the atoms 
fly apart. In harmonically trapped Bose-Einstein con- 
densates, the cloud starts to rotate when fl = ilc, and 
a vortex state forms at its center. As fl increases, more 
vortices appear and eventually form an array. As indi- 
cated, fl cannot exceed the trap frequency, w. 

Theoretical studies suggest that this picture becomes 
richer when an anharmonic term is added to the trapping 
potential [6-14]. The presence of the anharmonic term 
leads to two significant difi^erences. First, r2 is no longer 
bounded by u>. Second, the system can exhibit multiply- 
quantized vortex states, i.e., vortex states with more than 
one units of circulation (as opposed to singly-quantized 
vortex states which have one unit of circulation.) In- 
deed, for fixed interaction strength and sufficiently large 
n, the gas always exists in a state of multiple quantiza- 
tion. In the opposite extreme of fixed f2 and large interac- 
tion strength, the gas forms an array of vortices as in the 
case of harmonic confinement. In the intermediate regime 
where both parameters are large, there is a third phase 
which contains both an array of vortices and a strong sup- 
pression of the density in the middle of the cloud. This 
state is a mixed phase with a multiply-quantized vortex 
state at the center of the cloud and singly-quantized vor- 
tices around it. Recently, the experiment of Ref. [15] 
managed to create vortices in an anharmonic trapping 



potential and investigated the phases of the system as a 
function of il. Reference [16] has examined the dynamics 
of the gas in the presence of vortices. 

The purpose of the present study is to offer a descrip- 
tion of the transitions between phases of multiple and 
single quantization. These transitions can be of either 
first or second order. In the limit of weak coupling and 
small anharmonicity, this will lead to an exact descrip- 
tion of the continuous (i.e., second-order) transitions. In 
the same limit, the phase diagram as a function of the 
strength of the interaction versus will be seen to be 
universal. This phase diagram turns out to be very rich. 
Of particular importance are the triple points which we 
discover, which indicate the meeting of one phase of mul- 
tiple quantization with two other phases which are dis- 
tinct combinations of singly- and multiply- quantized vor- 
tex states with discrete rotational symmetries. In fact, 
we encounter an entire family of triple points. 

We present our model in Sec. II. We then examine 
separately continuous phase transitions in Sec. Ill and 
discontinuous phase transitions in Sec. IV. In Sec. V, 
we demonstrate the universality of the phase diagram 
and the location of the phase boundaries describing con- 
tinuous transitions in the limit of small anharmonicity 
and weak coupling. We offer a physical interpretation 
of our results in Sec. VI. Section VII contains the phase 
diagram that results from the solution of the exact eigen- 
value problem for the anharmonic confining potential and 
Sec. VIII discusses the direct numerical solutions of the 
Gross-Pitaevskii equation, which are in agreement with 
our variational approach. Finally Sec. IX gives a sum- 
mary of our results and conclusions. 



II. MODEL ~ APPROACH 

We start by assuming a confining potential of the form 

(1) 



F(p) = iMc.V[l + A(^)^] 
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Here p is the cylindrical polar coordinate, is the atomic 



mass, oo 



[h/Mu 



Nl/2 



is the oscillator length, and A is a 



small dimensionless constant. In the experiment of Ref. 
[15], A « 10~^. We neglect the trapping potential along 
the z axis since the cloud rotates about this direction 
and instead assume a constant density per unit length 
a = N/Z, where iY is the number of atoms and Z is the 
width of the gas along the z axis. The energy levels of 
the non-interacting (two-dimensional) problem are thus 
given by 



e = {2nr + \m\ + 1)?ujJ, 



(2) 



where rir is the radial quantum number, and m is the 
quantum number corresponding to the angular momen- 
tum. This equation emphasizes the large degeneracy 
which is associated with the many different ways of dis- 
tributing L units of angular momentum to A'' atoms. As 
shown in Refs. [17,18], this degeneracy is lifted by the in- 
teractions. One of the basic conclusions of these studies is 
that repulsive interactions always favor singly-quantized 
vortex states. On the other hand, when A > and the 
interaction is sufhciently weak, it is energetically favor- 
able for the system to form multiply-quantized vortex 
states. As the interaction increases, the effective repul- 
sion between vortices ultimately results in the splitting 
of multiply-quantized vortex states into singly- quantized 
states [6-14]. 

Here we shall account for the interactions with the 
usual assumption of a contact potential 



(3) 



where C/q = 47r?i a/M is the strength of the effective; two- 
body interaction with a equal to the scattering length 
for atom-atom collisions. The interaction is assumed 
repulsive, a > 0. (Reference [11] has examined an 
effectively-attractive Bose-Einstein condensate confined 
in an anharmonic potential.) The dimensionless quantity 
which plays the role of a coupling constant in our two- 
dimensional problem is thus aa, since the typical atom 
density n is ~ N/ {-Ka^Z) in the limit of weak interactions 
that we consider here. Therefore, the typical interaction 
energy uUq \s ^ aahu. 

As shown in Ref. [14], the phase transition between 
multiple and single quantization occurs when era is ~ 
A/a, where a is a dimensionless constant of order 10~^. 
Given that A w 10^"^ in Ref. [15], aa is less than unity in 
the region of the transition. This fact allows us to adopt 
a variational approach where the unperturbed states are 
restricted to the nodeless cigcnstates of the harmonic po- 
tential with angular momentum mfi 
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Here (j) is the angle in cylindrical polar coordinates. In 
this basis we can write the order parameter as 



*(p,0) = ^Cm^„ 



(5) 



where the coefficients Cm are variational parameters (as- 
sumed to be real without loss of generality). As a first 
step, we calculate the energy. From Ref. [14] we find that 
the energy per particle in the rotating frame in the state 
given by Eq. (5) is 

1 ^ , 
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J2 X] (^rnCnClCk{m,n\Vint\l,k)dm+n,l+k, 
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is the single-particle energy in the rotating frame mea- 
sured with respect to the zero-point energy, ?iu). Here the 
matrix elements of the interaction are given as 



{m,n\Vint\l,k) = aa 



(H + H)! 



2l"l+l»lv'|m|!|n|!|;|!|fc|! 



^m-\-n,l-\-k 



■ (8) 



As shown in Ref. [14], for sufficiently small aa this 
energy is minimized when only one component m = mo 
in the above sum is nonzero, i.e., Cm = for m 7^ mo. 
The critical frequencies, which denote the lower limit on 
the absolute stability of the state mo, are given by [14] 



mo 
\mo\ 



[l + A(lmol + l) 
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(2lmol -2)! 
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(9) 



The (straight) solid lines in Fig. 1 show the phase bound- 
aries in the fl/u - aa plane as given by Eq. (9). The 
transitions across these phase boundaries are discontinu- 
ous and of first order. 

As aa increases, however, the states with a single com- 
ponent m = mo become unstable resiilting in both con- 
tinuous as well as discontinuous phase transitions. In 
Sec. HI we identify the phase boundaries within a given 
sector mo as continuous transitions which are denoted as 
dashed lines in Fig. 1. Clearly, so long as aa is smaller 
than the value for which a given solid line is first cut 
by a dashed line, the straight solid lines correctly locate 
the discontinuous phase transitions between multiply- 
quantized vortices. On the other hand, for somewhat 
larger values of aa there is a competition between the 
mixed state that lies to the right of each solid line and 
the multiply-quantized state to its left. As a result of this 
competition, the straight lines of Eq. (9) no longer de- 
scribe the phase boundary correctly; the correct bound- 
ary is pushed to the left as shown in the dotted-dashed 
lines in Fig. 1. The determination of this boundary is 
given in Sec. IV. 
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FIG. 1. The phase diagram in the Q,/u! - aa plane for 
A = 0.05. The straight Mnes show the phase boundary be- 
tween multiply-quantized vortices given by Eq. (9) , which de- 
note discontinuous transitions. The dotted lines give the 
phase boundary between mo = 2 and (mi, mo, m2) = (0, 2, 4) 
(left), and between mo = 3 and (mi,mo,m2) = (0,3,6) 
(right), which denote continuous transitions. The dot- 
tcd-dashed lines show the phase boundary between mo = 1 
and (mi,mo,m2) = (0,2,4) (left), and between mo = 2 and 
(0,3,6) (right), which denote discontinuous transitions. The 
straight lines do not have any physical meaning for aa higher 
than the point where they are cut by the dashed lines, as the 
phase boundary is given by the dotted-dashed lines. 



III. CONTINUOUS TRANSITIONS 

We begin our consideration of the continuous transi- 
tions by showing that along the phase boundary, where 
only the m = mo state is present (and thus Cm — for 
m mo), the derivative of the energy with respect to 
all the Cm vanishes. Having established that this pure 
state is an cxtrcmum, wc will then examine the matrix 
that results by calculating the second derivatives of the 
energy with respect to any Cm and c„. The criterion for 
the stability of the pure state mo is then the positivity of 
the eigenvalues of this matrix. An instability and, hence, 
the phase boundary occurs when one eigenvalue of this 
matrix becomes negative. 

It is important to note that close to the phase bound- 
ary one needs to keep only those terms in the interaction 
energy [i.e., in the last term of Eq. (6)] which are at most 
bilinear in the Cm for m ^ mo- This implies that the 
only two terms which must be retained arc those with 
(mo, m,j|Vi„t|mo,m,,) (z = 1,2) and (mo, mo|l^nt|TOi, TO2) 
with mi + m2 = 2m,o. The second term is the only off- 
diagonal element in the second- derivative matrix. This 
observation simplifies the problem significantly since the 
matrix of second derivatives is block diagonal with the 




FIG. 2. The phase diagram of vortex states in a 
quadratic-plus-quartic potential in the Q/ui - aa plane for 
A — 0.05 with mo = 0,. . . ,15. The solid lines denote discon- 
tinuous transitions, and the dashed lines continuous transi- 
tions. The first circle on the left indicates the triple point, 
and the other circles show triple points on the lines of contin- 
uous transitions for mo = 8, 10, 13. 

dimensionality of each block being equal to two at most. 

The problem reduces to diagonalizing 2x2 matrices. 

Given the form of the interaction, Eq. (8), it is ele- 
mentary to verify that all derivatives dE/dCm do indeed 
vanish at the phase boundary as a consequence of an- 
gular momentum conservation [i.e., the delta function in 
Eq. (8)]. Turning to the second derivatives, d'^E/dcmdcru 
explicit calculation reveals that d'^E/dc^^ = and that 

-—^ = 2?iw{emi - emo) 
-|-[8(mo,mi|Mnt|TOo,'Ti,) - 4(mo, mo|yi„t|mo, mo)], (10) 
and 

= 4:{mo,mo\Vint\mi,m2)d2mo,mi+m2- (11) 

Using Eqs. (8), (10), and (11), we can now construct and 
diagonalize the resulting matrix, which is block diagonal 

as a consequence of the delta function in Eq. (11). Ev- 
idently, the highest dimensionality of each block is 2 x 
2. The result is shown as the dashed curves in Fig. 2 
for mo < 15. For mo < 6, we find that the most 
unstable mode is that which involves the states with 
{mi, mo, 1712) = (0,mo,2mo). However, as and mo 
increase, it becomes energetically unfavorable to have a 
significant atom density close to the center of the cloud. 
This is due to the form of the effective potential 
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FIG. 3. The same graph as in Fig. 2 axound the phase 
boundary for mo = 8. Sohd lines denote discontinuous tran- 
sitions and dashed lines continuous. The slope of the dashed 
curve changes at the center of the circle, which is a triple 
point. In the small triangular region the state is dominated 
by the (mi, mo, m2) = (2, 8, 14) components, while above the 
circle it is dominated by the (1,8,15) components. On the 
left the states arc mo = 7 (bottom), and (1,7, 13) (top). On 
the right the states are mo = 9 (bottom), and (2, 9, 16) (top). 

VM = lM{oj'-n')p' + ^Mcv^-^, (12) 

which has the famihar "Mexican hat" shape for O/w > 1. 
Since $o is the only state which does not vanish at 
the origin, this imphes that the admixtures of $o to 
the ground state wave function will become less favor- 
able as fi/w increases. Ultimately, the leading instability 
will involve states with mi ^ 0, and the ground state 
will have a hole at its center with a density that van- 
ishes like p^™!. Our results show that this happens at 
mo = 7 where the most unstable mode involves the states 
(mi,mo,TO2) = (1,7,13) instead of (0, 7, 14). These 
states lead to the most unstable mode over the entire 
region defined by Eq. (9). 

In the region where the toq = 8 mode is stable, the pic- 
ture becomes even more interesting. Initially, the most 
unstable mode involves the states with (TOi,mo,m2) = 
(2,8,14). When fi/w = 1.4240 and aa = 2.6880, how- 
ever, the state with (1, 8, 15) becomes the most unstable. 
This implies that there is a genuine triple point involv- 
ing the pure multiply-quantized vortex state with mo = 8 
and two phases which both have a vanishing density at 
the origin but which have different symmetries. As we 
shall see, the three phases that meet at the triple point 
correspond to a single multiply-quantized vortex, a lat- 
tice of singly quantized vortices, and a multiply quan- 
tized vortex surrounded by singly quantized ones. Fig- 
ure 3 shows these results. The circle denotes this genuine 



triple point. Note that the slope of the phase boundary 
(i.e., the dashed line) changes at this point. 

For higher values of mo this pattern persists with tran- 
sitions from some toq to either one set of (TOi,rno,m2) 
extending over the whole region of stability of the 
rrio state or with a transition from toq to one set of 
states (mi,mo,m2) initially, and then to some other 
set {m[,mf),m'2). The second alternative is realized for 
mo — 10 [which is initially unstable to (3, 10, 17) and fi- 
nally to (2, 10, 18)] and for mo = 13 [which is initially 
unstable to (5,13,21) and finally to (4,13,22)]. These 
triple points are indicated by the final two circles in 
Fig. 2. For mo < 15, these are the only two patterns 
found. Note that we cannot exclude the possibility of 
more triple points along a given segment of the second- 
order line for larger values of mo- In addition, mi is found 
to be an increasing function of mo. We off'er a physical 
interpretation of these results in Sec. VI. 



IV. DISCONTINUOUS TRANSITIONS 

Turning to the discontinuous transitions, it is instruc- 
tive to consider a specific example between the sectors 
with mo = 1 and mo = 2 (see Fig. 1). The value of aa at 
which the first dashed line cuts the solid line is 0.4226. 
For aa < 0.4226 the solid line thus gives the correct 
phase boundary between the pure states with mo = 1 
and m,o = 2, and the transition is discontinuous. For 
aa > 0.4226, there is a competition between the pure 
state = $1 and the state = co$o + C2$2 + C4$4 iden- 
tified in Sec. Ill (to the left of the dashed line) . Compar- 
ison of the energy of the system in the rotating frame in 
these two states allows us to find the boundary between 
the two phases (dotted-dashed line). Since each state de- 
scribes a local energy minimum, the transition between 
them is still discontinuous and of first order. Beyond 
aa — 0.4226 the solid line is thus meaningless, and the 
dotted-dashed line determines the phase boundary. Note 
that the slop(^ of tlw^ solid line is constant while that of 
the dotted-dashed line is not. This implies a discontinu- 
ity in the curvature at the joining of these curves. The 
use of a better variational wave function on the right of 
the solid curve indicates that the correct phase boundary 
must be pushed to the left. This is clearly seen in the 
results of Figs. 1, 2, 3, and 5. 

As mentioned above, the triple point, i.e., the point 
where the three phases coexist [10], lies on the bound- 
ary between the states with mo = 6, {mi, mo, 1712) = 
(0, 6, 12), and (1, 7, 13). To determine its precise location 
we need to combine the calculations presented in Sees. Ill 
and IV. This triple point is found to be at fi/w = 1.3633 
and aa = 2.1300 and is indicated by the first circle in the 
phase diagram of Fig. 2. 
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V. UNIVERSALITY AND EXACTNESS OF THE 
PHASE BOUNDARY 

The phase diagram of Fig. 2 is universal in the limit 
of weak coupling and small anharmonicity. More pre- 
cisely, if given values of aa and 1 — ^/u locate a point 
on the phase boundary of Fig. 2 for a specific A, the cor- 
responding values become {aa)' = f3aa and (1 — il/Lo)' = 
/3(1 — il/cj) for another A' = 6X. As A changes, both axes 
scale by the same amount and in that sense the phase di- 
agram of Fig. 2 is universal. Of course, this conclusion 
is valid only if the anharmonicity and the interaction are 
both sufficiently weak to justify our restriction to radial 
states with rir = 0. 

The method adopted in the present study of the con- 
tinuous transitions is also exact in the limit of weak 
coupling and small anharmonicity since, infinitesimally 
close to the corresponding boundaries, one can neglect 
all terms with m ^ mg and m 7^ in the expansion 
of Eq. (5). (In this regard, see also Sec. VI.) As one pro- 
ceeds away from the phase boundary, additional angu- 
lar momentum states must of course be included. With 
increasing coupling but arbitrarily close to the phase 
boundary, it may be necessary to consider states with 
rir 0. However, the restrictions on the angular mo- 
mentum can be maintained. Thus, our approach can 
be used to explore the phase boundary for significantly 
larger values of O/w and aa by including additional ra- 
dially excited states as needed. The phase boundary be- 
tween multiply-quantized vortex states is also calculated 
exactly in the same limit close to the line of second- 
order transitions. On the other hand, the calculation 
of phase boundaries for the discontinuous transitions in- 
volving multiply- quantized vortices and the mixed phases 
(i.e., the case considered in Sec. IV) is only approximate 
since additional angular momentum states are in princi- 
ple required for the description of the mixed phase. 

VI. PHYSICAL INTERPRETATION OF THE 
RESULTS 

It is useful to examine the spatial distribution of vor- 
tices and their multiplicity. The order parameter in the 
phases of multiple quantization is ^' = ^mo-, and just 
above the phase boundaries it has the form 

V]/ _ C„ /V I ^ fj JTOo I ^ fj xmi\ -\i\^ jlal 

(13) 

where mi < too < "^2 and z = pe'* = x+iy. A numerical 
calculation reveals that close to the phase boundary and 
to leading order, 

Cmi « Crn2 OC (5^^^, - 1 a 6^^'^, (14) 

where 6 = aa — {aa)c- Here {aa)c is the value of aa on 
the phase boundary for a fixed O. It is straightforward 



FIG. 4. Contour plots of the density of the cloud. The 
width in both axes is from — 4ao up to 4ao. The white color 
denotes high density and the black low density. The upper 
figures show the density of the cloud in the state with only 
mo = 6 (left) and (mi,mo,m2) — (0,6, 12) (right). The mid- 
dle figures refer to the case with only mo — 7 (left) and 
(1,7,13) (right). The bottom figures refer to the case with 
only mo = 8 (left), (2,8,14) (middle), and (1,8,15) (right). 
For mo > 7 there is always a node in the density at the center 
of the cloud. 

to develop a power-series expansion in S for all coeffi- 
cients that are not strictly zero (as a consequence of an- 
gular momentum conservation) [18]. For example, the 
next coefficient Cm contributing to the state of Eq. (13) 
is determined by the condition c^^^c^ ~ CmiCmoCm2Cm 
where m = mo -|- m2 — mi. This implies that 

(15) 

The next coefficient Cm' is determined by c^^c^, 
Crn„Cm2CmCm' , whcrc m' = m2 "I" m - mo = 2m2 - mi, 
and thus 

Cm' OC Cm^Cm OC <5^/^. (16) 

Evidently, the higher coefficients can be neglected at the 
phase boundary. 

Returning to Eq. (13), ^' has m2 nodes in the xy plane, 
and each of them corresponds to a vortex state. How- 
ever, only the mo of them are located in a region of non- 
negligible density since the density drops rapidly away 
from the center of the cloud due to the exponential fac- 
tor. More precisely, because of the first term on the right 
of Eq. (13), at a small 5 above the phase boundary there 
is still a multiply-quantized vortex state at the center of 
the cloud with mi units of circulation. In addition, there 
are mo — mi simple zeros or singly-quantized vortices 
surrounding it at a distance Ri which is 

:^ ^ f^V^''^""'"™'' (17) 
ao V "ii! / 
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The remaining vortices are located at a distance i?2 



R2 
ao 



moi 



l/[2(rn2-mo)] 



(18) 



For example, when the TOq = 8 state becomes unstable 
to mi = 2, mo = 8, and m2 = 14, there is a doubly- 
quantized vortex at the center of the gas, six more vor- 
tices around it at a distance of order Ri/ao ^ 2.35^^^"^ 
and a final six at a distance of order R^/ao ~ 3.4/5^/^^, 
where the density is essentially zero for sufficiently small 
5 (sec the bottom graphs in Fig. 4). Since m,i increases 
with increasing mo, just above the phase boundary a 
multiply-quantized vortex state splits into a multiply- 
quantized vortex state with lower circulation, with the re- 
maining (singly-quantized) vortices surrounding it. This 
is a mixed phase, and the hole which develops in the mid- 
dle of the cloud is a multiply-quantized vortex state with 
a quantum number mi, in agreement with Refs. [6,8]. 
It is worth mentioning that the maximum of the func- 
tion occurs for p/ao = \/m + 1, which gives the 
typical size of a vortex state with m quanta. 

Figure 4 shows contour plots of the density of the cloud 
for mo = 6, 7, and 8, and summarizes the results men- 
tioned above. These graphs all have distinct discrete ro- 
tational symmetries. Obviously, there cannot be a con- 
tinuous transition from one discrete symmetry to another 
unless one passes through the cylindrically symmetric 
triple point. Thus, such transitions are in general of first 
order. 
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FIG. 5. The same graph as in Fig. 2 for < mo < 8, with 
the single-particle eigenvalue problem solved exactly. The 
circle denotes the triple point. 



between the phases with mo = 7, (mi,mo,m2) = 
(0,7,14), and (1,8,15), and is shown inside the circle 
in Fig. 5. 



VIII. NUMERICAL SOLUTIONS OF THE 
GROSS-PITAEVSKII EQUATION 



VII. EXACT TREATMENT OF THE 
SINGLE-PARTICLE ENERGIES 

While our method treats the effect of the anharmonic 

term in the trapping potential perturbatively, it can eas- 
ily be extended to higher values of A and/or mo by em- 
ploying exact single-particle energies and wave functions 
which include the effects of the anharmonicity to all or- 
ders. 

In Fig. 5 we show the phase diagram that results from 
the use of such exact solutions for the non-interacting 
system (i.e., by solving the Schrodinger equation in the 
trapping potential V{p) rather than treating the anhar- 
monic term perturbatively.) The results are shown for 
A = 0.05 and mo < 8. In this case the phase bound- 
ary has no simple scaling properties and depends on the 
specific choice of this parameter. Although it is tempt- 
ing to regard this value of A as small, we find signifi- 
cant deviations from the results of Fig. 2. [Note, how- 
ever, the quadratic dependence of the perturbative en- 
ergy on m in Eq. (7) , which implies that perturbation 
theory becomes less reliable with increasing m]. In this 
case the most unstable mode for mo < 7 is of the form 
(mi,mo,m2) = (0, m,o,2mo), while for mo = 8 it is the 
one with (1,8, 15). Therefore the triple point occurs 



In addition to the above variational calculations, we 
have also solved the Gross-Pitaevskii equation numeri- 
cally. In order to locate the continuous phase transition 
lines to a few digits precision, we have used the follow- 
ing procedure. For given values of ua and fi/w, we take 
as initial state the harmonic-oscillator eigenfunction 
with the addition of small-amplitude random noise in the 
core region. (The appropriate value of mo is easily de- 
termined numerically provided that we are not close to 
the first-order transition and is consistent with the vari- 
ational results for the values reported here.) This wave 
function is propagated in imaginary time subject to the 
Gross-Pitaevskii equation in a rotating frame, and the 
overlaps, Cm', of the wave function with the harmonic- 
oscillator functions ^'m' are monitored for a range of m! . 
A steady increase in |cm'| for some m' indicates that 
the multiply-quantized vortex state is energetically un- 
stable towards a splitting into several vortices. If, on the 
other hand, all |cm'| decrease with time, the multiply- 
quantized vortex is the energetically favorable configura- 
tion. This procedure makes it easy to locate the position 
of the second-order transition without actually having to 
compute the ground state for each point in phase space. 
The latter procedure, while straightforward, is ill-suited 
to the problem at hand because convergence towards the 
ground state is extremely slow for the weak couplings 
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considered here. In adition, the difference in energy be- 
tween the many local energy minima is extremely small 
and requires excecdingiy high precision. In contrast, the 
method we have used is a fast and reliable way to lo- 
cate the continuous phase transition lines and represents 
a significant extension of the variational study. 

The results of this numerical study are shown in Fig. 6. 
The anharmonicity is chosen as A = 0.005. According to 
the universality shown in Sec. V, wc need only rescale 
the axes of Fig. 2 in order to obtain the variational phase 
diagram for this new value of A. Agreement with the vari- 
ational study is good for this very weak anharmonicity, 
but the discrepancy grows rapidly with increasing A. Ex- 
act wave functions have been computed for a few values 
of a a and and their symmetries are in agreement 

with those shown in Fig. 4. 

The first triple point occurs between mo — 6 and nio = 
7, precisely as predicted by the variational method in Sec. 
III. The most unstable mode for mo = 7 thus involves 
the states with (1,7,13). This is not inconsistent with 
the results reported in Sec. VII computed using exact 
anharmonic-oscillator eigenfunctions, because the value 
of A was different in the two cases. On the other hand, 
the present anharmonicity A = 0.005 is still large enough 
to result in one important difference compared to Figs. 
2-4. The triple point on the line of continuous transitions 
for Too = 8 is absent. Thus, the most unstable mode is 
the one with (toi, too, TO2) = (1, 8, 15) for all values of ft 
considered. We find instabilities of the form (toi, toq, TO2) 
with TOi > 2 only for too > 9. This again emphasizes 
that very weak anharmonicitics arc required to ensure 
the quantitative validity of the variational approach. 

IX. SUMMARY 

In summary, we have studied a rotating Bosc-Einstcin 
condensate with repulsive forces that is trapped in a 
quadratic-plus-quartic potential. As the rotational fre- 
quency of the cloud and the coupling between the atoms 
vary, the system exhibits three phases: a phase of 
multiple-quantization, a phase of singly-quantized vor- 
tices, and a mixed phase. Our calculated phase diagram 
turns out to be universal and partially exact in the limit 
of weak coupling and small anharmonicity. 

The authors wish to thank Gordon Baym, Ben Mottel- 
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dation. 
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FIG. 6. Same graph as in Fig. 2, for A = 0.005. The lines 
are identical to those in Fig. 2, except that the axes have 
been scaled according to the universality of the phase dia- 
gram. The crosses represent points of the continuous transi- 
tions computed numerically from the Gross-Pitaevskii equa- 
tion (see Sec. VIII). 
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